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Abstract 

Large-scale structures, observed today, are generally believed to have grown from ran- 
dom, small-amplitude inhomogeneities, present in the early Universe. We investigate 
how gravitational instability drives the distribution of these fluctuations away from the 
initial state, assumed to be Gaussian. Using second order perturbation theory, we calcu- 
late the skewness factor, S3 = (<5 3 ) / (<5 2 ) 2 . Here the brackets, (. . .), denote an ensemble 
average, and 5 is the density contrast field, smoothed with a low pass spatial filter. We 
show that S3 decreases with the slope of the fluctuation power spectrum; it depends 
only weakly on f2, the cosmological density parameter. We compare perturbative cal- 
culations with N-body experiments and find excellent agreement over a wide dynamic 
range. If galaxies trace the mass, measurements of S3 can be used to distinguish models 
with Gaussian initial conditions from their non-Gaussian alternatives. 
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1. Introduction 

Long before first estimates of the skewness in counts of galaxies became avail- 
able, Peebles (1980) showed how gravity can generate skewness in a random, initially 
Gaussian-distributed density field in an Einstein-de Sitter cosmological model. We have 
recently extended this calculation to Q / 1 (Bouchet et al. 1992, hereafter Paper I). 
In this Letter, we develop the formalism further in order to bridge the gap between 
theoretical concepts and observable quantities. We distinguish the mass density con- 
trast field, Sp/p, from the spatially smoothed field, S. While the former is not directly 
observable, the latter may be. Indeed, if the galaxies trace the mass, the moments of 
the frequency distribution in the counts of galaxies provide weighted averages of £ m - 
the m-point density correlation functions*, 

(n= f ^ 1 ;: (l% ux 1 ,...,x m ), a) 



rather than the proper moments, ((Sp/p) m ) = £ m (0). Here Xj are comoving spatial 
coordinates, dvi = F(-x.i)d 3 Xi, while F is a filter that determines the shape and volume, 
v = J d 3 xF(x), of a resolution element. In order to take the smoothing process into 
account, we calculate moments of a weighted average of Sp(x)/p, 

<5(x) = J ^(x') F(x - x') d 3 x' . (2) 

The filters considered here are spherically symmetric, sweep a unit volume and have a 
finite effective comoving half-width, R: 



I 



and jFl*)**,-*. (3) 



We derive analytic expressions for the skewness of the smoothed field 5 in models with 
various spectra of primeval fluctuations. We consider two kinds of filters, used by 
the observers and N-body simulators - the Gaussian and the "top hat", or a sphere 
with a fixed comoving radius. We show that the resulting skewness is sensitive to the 
slope of the spectrum of primeval fluctuations, as well as to the properties of the filter, 
contrary to incorrect claims made in the literature recently. We then test the limits of 
our approach and compare perturbative results with N-body experiments. We end by 
listing some still unresolved problems and summarizing our results. 



We neglect the shot noise terms. All our calculations are made in the fluid limit. 
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2. Perturbation theory 

Our calculations are based on the standard perturbative expansion for Sp/p in a 
Friedman universe, filled with non-relativistic pressureless fluid and zero cosmological 
constant (cf. §18 in Peebles 1980). We assume that to first order in perturbation theory, 
5p/p is a random Gaussian field. To first order, all its statistical properties are therefore 
determined by the power spectrum, P(k) = J £2(2) exp(ik • x) d 3 k/(2n) 3 . Perturbative 
calculations also assume that the amplitude of the density fluctuations, measured by 
their variance, is small: (S 2 ) <C 1. The linear order term in the perturbative expansion 
for (<5 2 ) is 




where Wk = f F(x) exp(ik • x) d 3 x is the smoothing window function. We consider 
three windows: Wk = 1 (no filtering at all), Wk = (3/kR) j\ (kR) , a top hat, and 
Wk = exp(—k 2 R 2 /2) , a Gaussian. Here and below ji denotes a spherical Bessel function 
of order I. Deviations from Gaussian behaviour, induced by gravity, appear at the 
second and higher orders and are fully determined by P(k). We define the skewness 
factor as the ratio of skewness to variance squared, S3 = (5 3 ) {d~ 2 )~ 2 ■ To lowest non- 
vanishing order, S3 is given by (Juszkiewicz & Bouchet 1991) 

/d 3 k d 3 V 
-^-^ P(k) P(k') W k W^ WVk'l T(K k') + 0{a 2 ) , (5) 

and the function T(k, k') is given by 

T(k, k') = 4 + 4/e(«) - 6 / u(fe/fc / ) + [2 - P 2 {p) , (6) 

where p, = k • k'/fcfc' , and P2 is a Legendre polynomial, while k is a slowly varying 
function of the current value of Q. For densities in the range 0.05 < SI < 3, 

n(n) w (3/14) n~ 2 / 63 , (7) 

(Paper I). The expression for (S 3 ) in k space, with fi = 1, but without smoothing, was 
first derived by Fry (1984), in a truly seminal paper. Goroff et al. (1986), whose analysis 
was also restricted to the Q = 1 case, were first to include the filters Wk- In the absence 
of filtering (Wk = 1), the dipole and quadrupole terms in equation (6) integrate to zero, 
and we obtain 

S3 = 4 + 4K(tt)«f + 6(^-2/63 _!) , ( 8 ) 

where the approximate form applies for 0.05 < Q < 3. The first term on the right 
hand side, 34/7, reproduces the the Peebles (1980) result for f2 = 1. The second term, 
found by Bouchet et al. (1992) is the "curvature correction", which arises when 0^1. 
This term is always small, and S3 is essentially insensitive to SI (this was pointed out 
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independently by Martel & Freudling 1991). To second order, the growth rates of 
the variance and skewness are such that for any comoving smoothing scale, the ratio 
(<5 3 ) / (S 2 ) 2 remains constant, if 0, = 1. In an open universe, S3 grows extremely slowly 
as Q decreases with time. In the SI — ► limit, k — > 1/4 (Paper I), and S3 — ► 5, a tiny 
increase compared to 34/7 pa 4.9. 

Equation (8) can be regarded as an approximation, valid in the regime when the 
density gradients across the filtering scale R are small. Indeed, let us consider a density 
field with a large coherence length, R c = o/o\, where 

a\ = j k 2 P{k) Wl d 3 k /(2Trf 

is the variance of the density gradient. In the limit R <C R c , we can Taylor ex- 
pand Wk W|k_k/| about the origin, k = k' = 0, and then evaluate the integral 
(5) termwise. Using the normalization conditions (3), and integrating over (i, we obtain 

S 3 = 4 + 4/e(0) - 2(R/R C ) 2 + 0{R/R c f . (9) 

The first two terms above describe the "local field" contribution, as in eq.(8). The third 
term is the tidal correction. Tidal effects, associated with the density gradients, tend 
to lower S3. This decrease is stronger for fields with smaller coherence length. In the 
next section we will see that this effect is also present for pure power-law spectra: S3 is 
anticorrelated with the relative amount of small-scale power. 



3. Power-law clustering models 

We will now study power law spectra, P oc k n . We set the smoothing length 
to unity and use dimensionless wavenumbers. Two kinds of filters are considered. We 
begin with the top hat. To separate the variables k and k' in the multiple integral (5), 
we use the addition theorem for Bessel functions *, 

jiM = £7 E^ + ^^^aT ' (10) 

where to = |k' — k|, and Pe((J,) are Legendre polynomials of order I. After substituting 
the above expansion to equation (5) and integrating over fx, the skewness factor can be 
expressed as 

S 8 = 4 + 4*(n)+f> + 40 (ii) 

* The formula for j\ can be derived by differentiating a similar expression for jo, given in 
Abramowitz & Stegun (1964, eq. [10.1.45]); for the full derivation, see Watson (1944). 
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where /3(n,£) = B 2 (n,£) if £ is even, and (3 = —B(n + l,£) B(n - 1,£) otherwise. The 
factors B(n,£) are integrals of spherical Bessel functions, 



ji(z)j e+1 (z) vrm!(^f±l)! 

For — 3 < n < 1, the series converges. Summing it requires some knowledge of the 
properties of the gamma function and a high threshold against boredom. The reward 
is the simple final result, 

S 3 = 4 + 4/e(«) - (3 + n) . (13) 

This expression is extremely weakly sensitive to U: the first two terms above are, in 
fact, identical to those on the right-hand side of equation (8). The skewness parameter 
is, however, strongly sensitive to the spectral slope, n. 

TABLE 1: Scale- free initial conditions 



B{- 



-m,£) = / 



n 


F(x) a 


S3 (perturbation theory) 


5 3 (N-body) 


a 


-3 


G 


34/7 = 4.9 






-2 


G 


(4vr + 9\/3)/7 =4.0 


3.8±0.2 e 


0.7 


-1 


G 


(20tt- 12>/3)/7V3 =3.5 


3.6±0.2 e 


0.7 





G 


8(4vrV3- 17)/7\/3 =3.1 


3.0±0.2 e 


0.7 


1 


G 


8(ll7r-12>/3)/2lV3 =3.0 






-1 


T 


(34/7) - 2 =2.9 


2.9±0.1 fe 


0.6 





T 


(34/7) - 3 =1.9 


1.8±0.3 C 


0.7 


1 


T 


1.9 d 


1.9±0.1 6 


0.7 



a G = Gaussian filter; T = top hat. 

b White (1992). 

c Bouchet & Hernquist (1992). 



Corrected for finite grid efects. 
e Weinberg & Cole (1992). 



We now consider the Gaussian filter, Wk = exp(— fc 2 /2). The variance and skewness are 
finite for n > —3. To evaluate the skewness integral (5) in this case, we used a coordinate 
transformation which reduces (k — k') 2 to a sum of squares. Similar multiple integrals 
are considered by Rice (1954). Our results for = 1 and several integer values of the 
power index n are listed in Table 1. In Figure 1 we plot 63(71) for Q = 1 and the entire 
range — 3 < n < 1, obtained by computing the integral (5) numerically. The case n = 1 
was considered earlier by Grinstein & Wise (1984). Again, when O^l, there is a small 
additional "curvature term" in the expression for S3. We write this additional term as 

AS 3 (n,n)^l (fr 2 / 63 -l)p„, (14) 
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so that the Gaussian-smoothed skewness factor is S^n, fi) = S^ra, 1) + AS3. For 
n = —3,-2,-1, and we obtained p n = 1, 1.04, 1.13, and 1.29, respectively. For 
n = —3, both (<5 3 ) and (<5 2 ) 2 diverge. However, their ratio is finite and equal to 4 + 4k. 
Note that in all cases, discussed in this section, the dependence of 53 on the relative 
amount of small-scale power is qualitatively similar to that in equation (9): S3 is a 
decreasing function of n. 



Figure 1. The skewness factor for power law spectra and Gaussian smooth- 
ing. Rigorous perturbation theory agrees well with the N-body results, while the 
Zel'dovich approximation systematically underestimates S3. 



4. Comparison with N-body experiments 

In Table 1 we summarize the results of several N-body experiments, used to study 
the evolution of skewness in the weakly nonlinear regime for scale-invariant initial con- 
ditions. The first column gives the power index, n. The second column specifies the 
spatial filter used in the simulation. The third column contains our predictions for S3, 
based on the perturbation theory. The fourth column gives a, measured at the same 
time and on the same smoothing scale in the simulation, as S3. The examples n = +1 
and —1 at the bottom of Table 1 were provided to us by Simon White, who calculated 
5*3, using the output from N-body simulations of Efstathiou et al. (1988). Similar results 
for 63 were recently obtained by Bouchet & Hernquist (1992) and Lahav et al. (1993) 
in their simulations of later stages of clustering (er « 1). For a Gaussian filter, we used 
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the recent numerical results of Weinberg & Cole (1992; we quote sampling errors as 
estimated by David Weinberg, private communication). 

We have also considered the cold dark matter model (regarded here as no more 
than an example of a spectrum with a scale-dependent slope). The skewness factor 
was calculated by numerical integration of eq. (5) and compared to an N-body simula- 
tion. For the simulation, we used the particle-in-cell code, described in Moutarde et al. 
(1991). This simulation involves 64 3 particles on a 128 3 grid. We assume the 'standard' 
spectrum with a Hubble constant of 50 km s _1 Mpc -1 , fi = 1, and no biasing. In 
Table 2 we list results for several smoothing widths R, and for two kinds of filters. All 
measurements were made at the time, when the r.m.s. fluctuation in the density field, 
smoothed with an R = 16 Mpc top hat, reached 0.92. 

TABLE 2: S 3 for a variable slope spectrum (CDM) 



R a 


F(x) b 


S3 (perturbative) 


S 3 (N-body) 


a(R) 


5.76 


G 


3.70 ± 0.07 c 


3.95±0.15 d 


1.30 


8.12 


G 


3.60 ± 0.07 


3.70 ± 0.20 


0.92 


11.52 


G 


3.50 ±0.07 


3.50 ±0.15 


0.64 


16.28 


G 


3.40 ± 0.07 


3.30 ±0.10 


0.44 


23.00 


G 


3.30 ± 0.07 


3.20 ±0.10 


0.29 


32.60 


G 


3.25 ±0.10 


3.05 ±0.10 


0.18 


11.52 


T 


3.15 ±0.10 


3.80 ±0.38 


1.30 


23.00 


T 


2.75 ±0.10 


3.05 ± 0.30 


0.64 


46.00 


T 


2.25 ±0.25 


2.45 ±0.25 


0.29 


65.20 


T 


2.00 ± 0.50 


2.20 ±0.20 


0.18 



a Smoothing scale in Mpc. 
b G — Gaussian filter; T — top hat. 
c Monte-Carlo integration errors. 
d Sampling errors. 

We were able to find only two cases of seemingly serious disagreement between pertur- 
bative predictions and numerical experiments. 

The first case is the n = 1 and a top hat filter, when the perturbative series (13) 
diverges: S3 = 00. Meanwhile, N-body experiments (Lahav et al. 1993, White 1992) 
give S3 « 2. This 'discrepancy' is easy to understand. The divergence is caused by 
the P(k) oc k behaviour at large wavenumbers. Such a spectrum is not reproduced 
in simulations for wavenumbers larger than the Nyquist frequency, fcjv, defined by the 
particle grid. To model this effect, we replaced the scale-free n = 1 spectrum, with 

ptjA _ J -^k) it k < kn] 
1 Akw, otherwise, 
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where A is a constant, and matches the grid used by Efstathiou et al. (1988). After 
this modification, the skewness integral (5) converged to S3 = 1.9, in perfect agreement 
with the N-body results (Table 1). 

The second case of seeming discrepancy is the S3 obtained by Park (1991) for 
n = — 1 and a Gaussian filter. The source of problem in this case is the Zel'dovich (1970) 
approximation (hereafter ZA), used by Park instead of an N-body code to calculate 
particle trajectories. ZA does not conserve momentum at second order in perturbation 
theory. This leads to an incorrect form for T(k, k'), with k in equation (6) set to 
zero (Grinstein & Wise 1986; Paper I). Using this incorrect expression, we were able 
to reproduce Park's result (Figure 1). Clearly, ZA systematically underestimates S3 
and disagrees rather badly with the rigorous perturbation theory and with the N-body 
results. 

5. Discussion 

Our results, summarized in equation (14) and Table 1, show that the skewness 
parameter is a decreasing function of the slope of the power spectrum; it is also sensitive 
to the shape of the window function used. 

Gravitationally induced skewness was also investigated by Coles & Frenk (1991, 
hereafter CF), who claim that S3 « 3 is a universal constant. The perturbative as well as 
N-body results, discussed here do not support this claim. The cause of our differences 
with CF is their failure to recognize that the perturbative results they use are not 
universal but instead are spectrum- and filter-specific. For example, their equation 
(11), quoted from Grinstein & Wise (1987, hereafter GW), is valid only for n = 1 and 
a Gaussian filter. More importantly, the GW formula is inapplicable to the N-body 
simulations, conducted by CF, because the spectra and filters assumed do not match 
each other; for the same reason the GW formula cannot be meaningfully compared 
with the S3, estimated from QDOT counts, considered by CF. Despite our quantitative 
differences, we do agree with CF on the qualitative level, that a scale-independent S3 
can be a signature of Gaussian initial conditions and a simple power-law spectrum (see 
the last paragraph in this section). 

We compared our calculations with results of N-body experiments to see if the 
perturbative series, truncated at second order, can lead to sensible results over a broad 
enough dynamic range. To our own surprise, the agreement between the two sets of 
results remained excellent even when the "small parameter" in the expansion, a, was 
close to unity. Apart from incorrectly conducted or misinterpreted simulations, all 
discrepancies appear to be within the sampling errors of numerical experiments. We 
also note that the agreement with N-body results was systematically better for the 
Gaussian filter than for the top hat, most likely because of high frequency sidelobes, 
which made the top hat-smoothed S3 sensitive to strongly non-linear fluctuations at 
small scales. 
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Remarkably, the qualitative properties of the weakly nonlinear clustering, de- 
scribed above, appear to hold in the nonlinear regime as well. When a S> 1, at least 
for = 1 and scale-free initial conditions, the relation between {S 3 } and {5 2 ) is well 
described by the semi-empirical formula 

£3 (xi,x 2 ,x 3 ) = Q [£ 2 (xi,x 2 )£ 2 (x 3 ,x 2 ) + sym.] , (16) 

where Q = Q(n) is a parameter, dependent on the initial spectral index only (Pee- 
bles 1980; Efstathiou et al. 1988). Adopting the Ansatz (16) is equivalent to setting 
T(k, k') = 3Q in the skewness integral (5). For a top hat filter, methods described in 
§3 yield 

a ~OQ 1 ^(3 + ^) 2 Q (m 

5s ~ 3g+ (5-^(2-^(3-^)4' (17) 

where v is the spectral index in the strongly-nonlinear regime, related to the initial 
slope by the so-called scaling solution, v = —6 / (n + 5) (Peebles 1980; if n > —3, then 
> v > — 3). Substituting Q(n), measured in N-body experiments (Efstathiou et al. 
1988), we get S3 = 4.5,2.9, and 2.3 for n = —1,0, and +1, respectively. Similar results 
were recently obtained in N-body experiments of Weinberg & Cole (1992), who directly 
measured Ss(n) in the nonlinear regime, and by Fry et al. (1993), who found that Q(n) 
decreases with n. To summarise: the scaling (5 3 ) oc u 4 seems to hold equally well both 
for (j < 1 and a S> 1; moreover, in both regimes S3 is a decreasing function of n. 

Preliminary results from observations appear to be consistent with a scale in- 
dependent S3 (Saunders et al. 1991; Coles & Frenk 1991; Park 1991; Bouchet et al. 
1991, 1993; Gaztahaga 1992; Lahav et al. 1993) This is exactly what is expected in the 
standard gravitational instability picture with Gaussian initial conditions and a sim- 
ple power-law spectrum. It is even possible that strongly non-Gaussian models, which 
give S3, diverging like cr _1 instead of being constant, can already be excluded (Silk 
& Juszkiewicz 1991). However, before reaching dramatic conclusions, more work is 
needed. For example, the absence of rich galaxy clusters in the IRAS survey is likely to 
cause a systematic underestimate of S3. We need to understand how S3 in the matter 
distribution relates to the skewness in galaxy counts, as matter and galaxies may be 
distributed differently. Finally, it is necessary to account for the effect of redshift space 
distortion on S3. Of all of the unresolved problems listed above, the latter admits the 
most straightforward solution, and we plan to report on this in near future. 

Acknowledgements: RJ thanks David Spergel and Chris McKee for raising his in- 
terest in deviations from Normal behaviour, and Alain Omont, John Bahcall and Simon 
White for their hospitality in Paris, Princeton and Cambridge, respectively. We are 
grateful to Ofer Lahav, David Weinberg and Simon White for sharing their unpublished 
numerical results with us. We also acknowledge partial financial support from the Pol- 
ish Council for Scientific Research (KBN), grant No. 2-1243-91-01. The computational 
means (CRAY-2) were made available thanks to the scientific council of the Centre de 
Calcul Vectoriel pour la Recherche. 



9 



References 



Abramowitz, M. & Stegun, I. A., 1964, Handbook of Mathematical Functions (Nationa Bureau 

of Standards: Washington). 
Bouchet, F.R., Davis, M, & Strauss, M., 1991, in "The Distribution of Matter in the Universe" , 

2nd DAEC Meeting, ed. G.A. Mamon & D. Gerbal (Meudon: Observatoire de Paris), 

p. 287. 

Bouchet, F.R. & Hernquist, L., 1992, Ap. J., 400, 25. 

Bouchet, F.R., Juszkiewicz, R., Colombi, S. & Pellat, R., 1992, Ap. J., 394, L5 (Paper I). 

Bouchet, F.R. et al, 1993, preprint. 

Coles, P. & Frenk, C, 1991, M.N.R.A.S., 253, 727. 

Efstathiou, G., Frenk, C.S., White, S.D.M. & Davis, M., 1988, M.N.R.A.S., 235, 715. 
Fry, J.N., 1984, Ap. J., 279, 499. 

Fry, J.N., Melott, A.L. & Shandarin, S.F., 1993, preprint. 
Gaztanaga, E., 1992, Ap. J., 398, L17. 

Goroff, M.H., Grinstein, B., Rey, S.-J. & Wise, M.B., 1986, Ap. J., 311, 6. 
Grinstein, B. & Wise, M.B., 1987, Ap. J., 320, 448. 

Juszkiewicz, R. & Bouchet, F.R., 1991, in "The Distribution of Matter in the Universe", 2nd 
DAEC Meeting, ed. G.A. Mamon & D. Gerbal (Meudon: Observatoire de Paris), p. 
301. 

Lahav, O., Itoh, M., Inagaki, S. & Suto, Y., 1993, Ap. J., 402, 387. 
Martel, H. & Freudling, W., 1991, Ap. J., 37, 11. 
Moutarde, F., et al. , 1991, Ap. J., 382, 377. 
Park, C, 1991, Ap. J., 382, L59. 

Peebles, P.J.E., 1980, "The Large Scale Structure of the Universe" '(Princeton: Princeton Uni- 
versity Press). 

Rice, S.O., 1954, in "Selected papers on noise and stochastic processes" , ed. N. Wax (New 

York: Dover), p. 73. 
Saunders, W., et al, 1991, Nature, 349, 32. 
Silk, J. & Juszkiewicz, R, 1991, Nature, 353, 386. 

Watson, G.N., 1944, A Treatise on the Theory of Bessel Functions (Cambridge: Cambridge 

University Press). 
Weinberg, D.H. & Cole, S., 1992, M.N.R.A.S., 259, 652. 
White, S.D.M. (1992), private communication. 
Zel'dovich, Ya.B., 1970, Astr. Ap., 5, 84. 



10 



